// ADDITIVE WAVETABLE SYNTHESIS
// D. TAYLOR 2014

#ifndef __WAVETABLE_H__
#define __WAVETABLE_H__

#include "config.h"
#include "intrinsics.h"
#include "fpmath.h"

// Wavetable interpolation method
#if WAVE_HQ_INTERPOLATE
#define WAVE_INTERPOLATE(A0, A1, A2, A3, R) (scubic((A0) << 14, (A1) << 15, (A2) << 15, (A3) << 14, R >> 1) << 4) // Cubic interpolation, 24bit result
#else
#define WAVE_INTERPOLATE(A0, A1, A2, A3, R) linear((A1) << 16, (A2) << 16, R) // Linear interpolation, 32bit result
#endif

// Readability macros
#define WAVE_SIZE  (1 << WAVE_BITS)         // Wavetable size

// Wavetable MIP-level calculation
__inline
uint32_t wave_calc_mip(uint32_t rate) {
  uint32_t mip = clz(rate);
  if(mip < 2) mip = 2;
  return mip > WAVE_BITS ? WAVE_BITS : mip;
}

// Wavetable MIP-level mu calculation
__inline
uint32_t wave_calc_mu(uint32_t rate) {
  uint32_t mip = clz(rate);
  if(mip <= 2 || mip > WAVE_BITS) return 0;
  return rate << (mip + 1);
}

// Sample wave table (with MIP-mapping)
__inline
samp32bit wave(uint32_t pos, uint32_t mip) {
  uint32_t register base = 1 << mip;
  uint32_t register index = pos >> (32 - mip);
  return wavetable[base + index] << 16;
}

// Sample wave table (with interpolation and MIP-mapping)
__inline
samp32bit fwave(uint32_t pos, uint32_t mip) {
  if(mip <= 2) return -(fastsin(pos) >> 1); // use sine for > nyquist
  uint32_t register base = 1 << mip;
  uint32_t register index = pos >> (32 - mip);
  return WAVE_INTERPOLATE(
    wavetable[base + ((index - 1) & (base - 1))],
    wavetable[base + index],
    wavetable[base + ((index + 1) & (base - 1))],
    wavetable[base + ((index + 2) & (base - 1))],
    pos << mip
  );
}

// Sample wave table (with interpolation and interpolated MIP-mapping)
samp32bit fmwave(uint32_t pos, uint32_t mip, uint32_t mu) {
  samp32bit sample = fwave(pos, mip);
  if(mip <= 2) return sample;
  return linear(sample, fwave(pos, mip - 1), mu);
}

// Halfband decimation filter, 5 taps
void halfband(samp16bit *out, samp16bit *in, uint32_t size) {
  int32_t z[5] = {0, 0, 0, 0, 0};
  int32_t x0, x1;
  for(int32_t n = -5; n < (int32_t)size; n++) {
    if(n < 0) {
      x0 = *(in + ((size + n) * 2) + 0);
      x1 = *(in + ((size + n) * 2) + 1);
    } else {
      x0 = *in++;
      x1 = *in++;
    }
    int32_t f[3] = {x0 * 852, x0 * -4167, x0 * 19699};
    if(n >= 0) out[(n + size - 3) & (size - 1)] = (z[4] + f[0]) >> 16;
    z[4] = z[3] + f[1];
    z[3] = z[2] + f[2];
    z[2] = z[1] + f[2] + (x1 << 15);
    z[1] = z[0] + f[1];
    z[0] = f[0];
  }
}

// Auto-generate MIP-maps by decimation
void build_mipmaps() {
  unsigned n;
  for(n = WAVE_BITS - 1; n > 1; n--)
    halfband(&wavetable[1 << n], &wavetable[2 << n], 1 << n);
}

#endif
